Mon. Not. R. Astron. Soc. 000, 000-000 (0000) 



Printed 2 February 2008 



(MN M?eX style file vl.4) 



Substructure and the halo model of large-scale structure 

Ravi K. Sheth^ & Bhuvnesh Jain^ 

^ Department of Physics and Astronomy, University of Pittsburgh, 3941 O'Hara Street, Pittsburgh, PA 15260 
^ Department of Physics and Astronomy, University of Pennsylvania, 209 S. 33 Street, Philadelphia, PA 19104 



2 February 2008 



(N 
O 
O 
(N 

< 
0^ 



> 
m 
in 
m 

00 

o 

(N 
O 

Oh 

6 



ABSTRACT 

We develop the formalism to include substructure in the halo model of clustering. Real 
halos are not likely to be perfectly smooth, but have substructure which has so far 
been neglected in the halo model — our formalism allows one to estimate the effects 
of this substructure on measures of clustering. We derive expressions for the two-point 
correlation function, the power-spectrum, the cross-correlation between galaxies and 
mass, as well as higher order clustering measures. Simple forms of the formulae are 
obtained for the limit in which the size of the substructure and mass fraction in it is 
small. Inclusion of substructure allows for a more accurate analysis of the statistical 
effects of gravitational lensing. It can also bring the halo model predictions into better 
agreement with the small-scale structure seen in recent high resolution simulations of 
hierarchical clustering. 
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1 INTRODUCTION 

Sheth & Jain (1997) described how the halo model for clus- 
tering can allow one to model the distribution of matter 
in the highly nonlinear regime. The model falls within the 
broader framework described by Neyman & Scott (1954) 
and Scherrer & Bertschinger (1991). It combines results from 
Peebles (1974) and McClelland & Silk (1977) with the work 
of Press & Schechter (1974), and is able to provide a good 
description of nonlinear clustering seen in numerical simula- 
tions of hierarchical gravitational clustering. The halo model 
assumes that most of the mass in the Universe is bound up in 
virialized dark matter halos, and that statistical measures 
of clustering on small scales are dominated by the inter- 
nal structure of the halos. The agreement with simulations 
shows that is possible to provide an accurate description of 
clustering in the small-scale nonlinear regime even if one has 
no knowledge of if and how the halos themselves are clus- 
tered. This halo-model of clustering has been the subject 
of much recent interest (e.g. Seljak 2000; Peacock & Smith 
2000; Ma & Fry 2000; Scoccimarro et al. 2001; Cooray & 
Sheth 2002). 

To date, almost all analytic work based on the halo- 
model approach assumes that halos are spherically symmet- 
ric, and that the density run around each halo center is 
smooth. Halos which form in numerical simulations of hi- 
erarchical clustering are neither spherically symmetric nor 
smooth (e.g., Navarro, Frenk & White 1997; Moore et al. 
1999; Jing & Suto 2002). About ten percent of the mass of 
a halo is associated with subclumps (Tormen, Diaferio & 



Syer 1998; Ghigna et al. 1999). The main purpose of the 
present work is to derive a model which accounts for this 
substructure. 

Section |^ shows how to compute two-point statistics, 
the correlation function and its Fourier transform, the power 
spectrum, when substructure is important. Section ^ shows 
that the model can be easily extended to estimate higher- 
order statistics. Section ^ shows how our formalism can in- 
corporate a range of parent halo and subclump masses, and 
Section |^ provides a few explicit examples. This section in- 
cludes a discussion of how to model the shape of the sub- 
clump mass function. Section ^ summarizes our results, and 
suggests various other applications of our formalism. 



2 TWO-POINT CORRELATIONS 

The halo model approach assumes that all mass is bound 
up in dark matter halos. All statistical measures of cluster- 
ing are then decomposed into two distinct types of contribu- 
tions: one comes from sets of particles which are in the same 
halo, and the other comes from particles which are in differ- 
ent halos. The assumption is that, on the smallest scales, it 
is the single-halo contribution which dominates the statis- 
tic. In the context of the present paper, this means that 
we expect that the contribution to, e.g., the power spec- 
trum, which comes from substructure cannot be important 
on scales larger than a typical halo. Therefore, even if sub- 
structure changes the single-halo contribution substantially 
compared to the case of smooth halos, it should be accurate 
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to ignore substructure when estimating the two- and higher- 
order halo terms which dominate on larger scales. The pri- 
mary consequence of substructure, then, is that the single- 
halo term becomes more complicated. For example, for two- 
point statistics, the one-halo term will now have three types 
of pairs: both particles from the smooth component, both 
particles from the substructure component, and one particle 
from the smooth component with the other from the sub- 
structure. Our goal will be to derive expressions for these 
different contributions. 



2.1 Preliminaries 

The probability of finding a subclump in the volume element 
(f'r at distance r from the center of the parent halo of mass 
M is 



J (Pr nc{r) 



(1) 



Here Uc describes the number density of subclumps at r. 
In what follows, it will prove convenient to multiply this 
number density by a mass and so define a mass density pc. 
We choose this mass so that J ct^r pc{r) — M, the total 
mass in the halo. (For simplicity we will sometimes assume 
that the subclumps have the same density run as the smooth 
dark matter distribution, but our analysis is not confined to 
this case. In this special case, pc is the same as the mass 
density profile of the smooth component.) 

The density at r from the halo center is the sum of 
the smoothly distributed mass plus the contribution from 
the subclumps. If F and fi denote the density run in units 
of the mean density p around the parent halo and the zth 
subclump respectively, then 



p{r) 



(2) 



The mass in the smooth and subclump components is de- 
fined by 



Ms 



drA-Kr^ p F{r) and 



dr4nr^ p f{r) 



(3) 



The total mass M is the sum of the smooth and clumped 
components. 

Our strategy will be to first derive expressions for, e.g., 
the correlation function for fixed values of M, rrii and A'', 
and to average over the distributions of these variables later. 
We will keep the discussion as general as possible; realistic 
choices for F , f, p, and the halo and subclump mass func- 
tions will be inserted into the formalism later. 



2.2 The correlation function of the mass 

The ensemble-averaged two-point correlation function ^(r) 
will be obtained by averaging over halos with distributions 
of sub-structure specified by p{r). We will first obtain the 
two-point correlation function within a halo, denoted C(r), 
then average over substrucutre to get S,{r), and finally in 
section 4 consider the averaging over varying numbers of 
subclumps and of parent halo masses. For a halo of total 
mass M we define 



C{r) 



d s 5[s)5{s + r), where 5{r) 



P 



1, (4) 



and n denotes the number density of the parent halos. If 
the mass densities within the subclumps are all much larger 
than the mean density p, then 5{r) ~ p{r)/p. Writing out 
all the terms explicitly shows that 



C(r) 



d^s F{s)F{s + r) 



+ d'^s fi{s - ri) fi{s + r - Vi) 

i 

+ ^ J d's F{s + r)f,{s^r,) 

i 

+ J2 J d'^s F{s)Ms + r-r,). (5) 

i 

This is the correlation function associated with a smooth 
component and the i subclumps at the specified positions. 
The first term is the contribution from particles which are 
not in the subclumps, and is familiar from previous work; 
this is the contribution from the first pair-type. The second 
term is from pairs where both particles are in the same sub- 
clump, whereas the third term is from pairs where the two 
particles are in separate subclumps. These two terms repre- 
sent the contribution from the second pair-type. The fourth 
and fifth terms are from pairs in which one particle is in a 
subclump and the other is in the smoother component. 

We are less interested in the correlation function associ- 
ated with a specific realization of the subclump distribution, 
than we are with what happens upon averaging over the var- 
ious possible subclump distributions. That is, we are more 
interested in 



d?r. 



d^n 



p{r,)p{rj) C{r) 



(6) 



The next step is to compute these averages. 

Let Ass(r) denote the value of the first integral in ex- 
pression It is the same for all the possible subclump 
distributions, so the result of averaging it over the distribu- 
tions is 



Xssir) = / d^sF{s)F{s + r) 



(7) 



This term is the convolution of F with itself, and is the only 
term which most halo-models use. The remaining terms are 
due to the substructure. Before we compute them, it will 
prove convenient to rewrite this first term as 




Here Asniooth(»') denotes the correlation function if the to- 
tal mass were smoothly distributed around the center of the 
parent halo (i.e, if there were no subclumps). Thus, the fac- 
tor of {Ms/MY simply denotes the fact that now only a 
fraction of the mass is in the smooth component. 
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Let Xu{r) denote the result of averaging each of the 
second terms in equation (m over p{ri). Then 



Ki{r) = / d?rip{ri) I cf s fi{s - ri) fi{s + r - ri) 



d X fi{x) fi{x + r), 



(8) 



where we have set x = s — ri. This term is the convolution 
of each subclump profile fi with itself. 

Averaging the third term in equation (^) is more com- 
plicated so we will do it last. The result of averaging the 
fourth term is 

Xs^{r) = d\,p{r,) d''sF{s + r),U{s-r^) 



dh I d\, Ms-r,) F(s + r) ^'^''''^ 



d^x f,{x) I dV, F{x + r,+r 



J d^r pc{r) 
' M 



d-^x 



M/p 



'- fi(x) Xsc{x + r) 



(9) 



where Xac denotes the convolution of F with Pc/p- (If the 
subclumps were distributed around the halo center similarly 
to the dark matter, i.e., Pc/p oc F, then Asc oc Ass.) By 
symmetry, the average of the fifth term is the same. 
And finally, the average of the third term is 



Aij(r) = I (fn I d:''rj p{ri)p{rj) 



X j d's f,{s~r,)f,{s + r~rj), (10) 

where the integral over s is the convolution of the two sub- 
clump profiles. Written this way, Aij is the weighted sum of 
this convolution over all pairs of subclump positions. It is 
interesting to re-arrange the order of the integrals above: 



Xij{r) = / d\i / d\jp{ri)p{rj) 



X d s f,Xs ~r,) f,{s + r - rj) 



d S d Ti fi{s — Ti) 



X I d rjfj(s - + r + rij)p{rij + rj)p{rj) 



d X d rij fi{x) fj{x + r + rij) 



[m/pY 



d'x j d?vH^)h{v) ^'\,/'ly (11) 

where Acc denotes the convolution of pc/p with itself. The 
second line follows from setting r^j = Ti — rj, the third line 
from setting a; = s — and the final expression from setting 
y = x + r + rij . When written in this way, Xij appears very 
similar to the two-halo term when substructure is absent; 
Acc plays the role of the 'subclump-subclump correlation 
function' (compare equation 2 in Sheth et al. 2001). 

Using the above results, the correlation function ^(r) 
defined in equation (^) can be expressed as 



= Ass(r)+^ A„(r)+2^ As.(r)+^ ^ A,,(r).(12) 

i i i j^i 

Note that all the A terms on the right-hand side have di- 
mensions of volume, so that ^ is dimensionless. 



2.3 Limiting cases for small-sized subclumps 

Now consider some simple limiting cases. If the subclumps 
are much smaller than the parent halo, then we can treat 
them as point masses when performing the integrals which 
define Asi and Xij. The associated delta function profiles 
simplify the integrals, so that 

C(r 



^ ^ Ass(r)+^A.,(r) + 2^^Asc(r) 

-EE^iM.). 



(13) 



i ii^j 



If in addition, the distribution of subclumps around the halo 
center is similar to the dark matter, pc oc F, then we can 
set Asc oc Ass and Acc oc Ass- In this case, the previous 
expression becomes 



mirrij 



i ji^i 



+ 



-^smooth 

(r) 



E 



(Recall that Asmooth denotes the correlation function if all 
the mass was smoothly distributed around the halo center.) 
If the fraction of the total mass which is in subclumps is 
smaU, then ^(r) ^ nAsmooth(r') + "nXair). In this limit, 
the total correlation function is well approximated by tak- 
ing what one would have got if the mass was smoothly dis- 
tributed, and then adding the contribution from the indi- 
vidual subclump components. In the small mass and size 
limit, the contribution from the subclumps can only be im- 
portant on scales smaller than the typical subclump, so that 
this additional contribution is only important on very small 
scales. On scales larger than the diameter of a typical sub- 
clump, the correlation function looks just as though the mass 
within halos is smoothly distributed. This simple assump- 
tion is probably sufficiently accurate for most applications. 

We remarked that the term Acc(t') could be thought of 
as the correlation function of the subclumps. To see why, 
suppose there is no mass in the smooth component, and all 
the subclumps are infinitesimally small: i.e., we replace all 
factors of the subclump density profile /i with delta func- 
tions &T>(ji)- Then, when r 7^ 0, only the third term in 
equation (^) contributes any pairs. Using the delta functions 
reduces this term to 



A».(r) = EE/ '^'^^ / 



d Tj p(ri)p(rj) 



i j^i 



X / d s ft{s - r,) fj{s + r - rj 
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EER 



d ri p{r,)p{ri + r), (15) 



which is proportional to the convolution of the subclump 
distribution pc with itself. 



2.4 The cross— correlation between subclumps and 
mass 

We can use a similar argument to compute the cross cor- 
relation between subclumps and mass. That is, we imagine 
we sit at the center of the ith subclump and we then com- 
pute the typical density at distance r from it. Because the 
first term in equation is not centered on a subclump it 
no longer contributes. Since our constraint only requires one 
member of each pair of positions to be centered on a sub- 
clump we now need to replace one factor of fi with a delta 
function. This means 



n Xii{r) 
n \ij{r) 



fi{r), n \, 



^cc{y - r) 



and 



(M/p)2 ' 

so that the subclump-mass cross correlation function is 



Cx(r) 



E 



M/p 



+ 



E/'(^ 



+ 



EE/ d'^fAs) 



Aee(5 - r) 

{M/py 



(16) 



In the limit in which the subclumps are much smaller in size 
than the parent halo the last term simplifies to give, 



M/p 

When pc/ p (X F this becomes 
Asa(r) 



M {M/p) 



M/p 



1+E 



M 



+ 



E/»( 



(17) 



and further, if the total mass in subclumps is small we obtain 
Cx (r) ^ ^ss{r)/{M/p) + Mr). (18) 

i i 

It is easy to verify that the result for can be obtained 
simply by averaging 

i i i 

+EE^^(^'+^"^^) (19) 

over all realizations of subclump positions, p{ri) and p{rj). 
2.5 Power spectra 

The power spectrum is the Fourier Transform of the cor- 
relation function. Since the correlation function involves a 
number of convolution-type integrals, the power spectrum is 
given by simple multiplications of the various density profile 
factors. Let 



U{k) 



J dr 47rr^ F{r) sin(fcr-) /kr 
(27r)3 J dr4-Kr'2 F{r) 



(20) 



denote the Fourier transform of the density run of the 
smooth component F normalized by the total mass con- 
tained in the profile, and define Ui{k) and Ucik) similarly (we 
have assumed we are working in three-dimensions). Then 

dr sin(fcr) 



P{k) 



(27r)3 


kr 




/Ms 




\ M 




+2 



'U{kf 
rrnMs 



■EE 



M2 

rriimj 
M2 



u,{k)U(k)Uc{k) 



Ui{k) Uj{k) Uc{k) 



+ 



(21) 



which we could write as 

P{k) = Pss{k) + 2 ^ P,,(fc) + ^ ^ P,,(fc) + J2 P"ik)- 

i i j^i i 

A similar calculation shows that the cross spectrum between 



subclumps and mass is 



— U{k) Uc(k) 
P 



+ (7,(fc)2^^^*,(fc)-t-^u.(fc) 
and the power spectrum of the subclumps is 

If we set Uc(k) = U{k), and further assume that the 
subclumps are much smaller than the parent halo, then the 
power spectrum reduces to the sum of the power spectra of 
the individual components: 



(22) 



(23) 



P{k) 



-Ely 



Ui{k) 



U{kf 



(24) 



This shows that whether or not the substructure component 
dominates the small scale power depends on the fraction of 
mass in the subclumps, and on how much more dense the 
subclumps are compared to the smooth component. 

To see this a little more clearly, consider a specific ex- 
ample. Suppose that the density runs are Gaussians, and 
that all the A'^ subclumps have the same mass rrii and size 
Ri. If i? denotes the characteristic scale of the density runs 
of the smooth component, then U = exp{—k^R^/2) and 
Ui = exp{—k^Ri/2). The subclumps dominate the power if 
N{mi/p)^Ui > {M/pY U{k)^. This happens on scales where 
fc2 > [21n(M/M,0+lnAr]/iiV[l-(-R«/-R)^]- where we have 
set Mci = Nrrii. Note that if the total mass in subclumps is 
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Figure 1. Effect on the power-spectrum as the parameters which describe the subclump distribution are varied. The dashed Hne, which 
is the same in both panels, shows A(fc) when all the mass is in smooth NFW halos of radius rvir within which the average density is 
200 times the background density. The solid curves in the panel on the left show A(A;) when / = 0.1 of the total mass of each parent 
halo is in subclumps, each containing 1% of the mass of the parent. The average density within each subclump is 500p (lower) and 5000p 
(higher). In the panel on the right, / = 0.9, the subclumps are 500 times denser than the background, and the subclumps are each 0.01 
(peak at larger k) and 0.1 (peak at smaller k) times the mass of their parents. 



small, Mci/M <C 1, then the subclumps dominate the power 
only at very large fc; if this mass is divided up among many 
subclumps, the subclumps only dominate at even higher k. 
For a fixed mass ratio, the exact scale on which the sub- 
clumps dominate depends on the size ratio Ri/R. This sug- 
gests that a feature in P{k) at large k may provide informa- 
tion about the nature of the subclumps. Simulations suggest 
that M/Mci ~ 10 and Ri ^ R/10, which yields a critical 
value of fc \/6 + lnN /R. It is interesting that this is just 
beyond the reach of current simulations. 

Fig. 1^ illustrates this with slightly more realistic Navarro, 
Frenk & White (1997) density profiles. The various curves 
show A(fc) = Ank^ P{k) for models in which the density 
field is made up of Poisson distributed NFW halos. The 
parent halos are truncated at their virial radii rvir defined 
so that the average density within Tvir is 200 times that 
of the background. We set the NFW core radius to be 
flNFW = 0.1 rvir- The dashed line (same in both panels) 
shows A(fc) if there is no substructure. We then assumed 
that a fraction / of the mass of each parent halo was in sub- 
clumps, each of mass m. For simplicity, we assumed that the 
distribution of subclumps around the halo center was given 
by the same NFW form, and that the distribution of parti- 



cles around each subclump center was also NFW, with core 
radius onfw = 0.1 rsc- The value of rac was set by requir- 
ing that the average density within the subclumps equal 5sc 
times the background densiU^ The solid curves in the panel 
on the left show equation ( |2l|) when / = 0.1, m = O.OIM 
and Ssc ~ 500 (lower) and 5sc ~ 5000 (upper). This shows 
that, all other things being equal, denser subclumps con- 
tribute more power. The solid curves in the panel on the 
right show results when f — 0.9, Ssc = 500 and m/M = 0.1 
and m/M = 0.01 (departure from dashed curve apparent at 
high and still higher k). Comparison with the panel on the 
left shows that increasing the fraction of mass in subclumps, 
/, increases the small scale power. The scale on which this in- 
crease becomes apparent depends on {m/M)^^''^{Svii/5sc)^^''^ , 
the typical radii of the subclumps. 



3 HIGHER-ORDER STATISTICS 

Higher order correlations at a given small scale are domi- 
nated more strongly by the one-halo term than is the two 
point function, so the effect of substructure on the one-halo 
term is of great interest for higher order correlations. In this 
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section we will consider the effect of substructure on the 3- 
point correlation function. The extension to higher orders is 
obvious. 

We begin with the expansion of the density p(r) given 
in equation The 3-point function in real space is then 
defined by 



^123 = {5{ri)5(r2)5{r3)) ~ — {p{r{)p{r2)p(r3)) 



(25) 



It is convenient to consider the expressions for higher or- 
der correlations in Fourier space, since we will find that in 
the limiting case of interest the only additional term due to 
substructure is a term involving products of the substruc- 
ture profile. 

The Fourier transform of the 3-point function is the 
bispectrum B defined by 



5{ki)5{k2)5{kz] 



B(fci,fe2,fe3) {2-kY Soiki 



(26) 



where the Delta function indicates that ki + k^ + ks = 0. 

It is the sum of contributions from triplets which are all 
in the smooth component, triplets all in the same halo, plus 
contributions from the various cross terms. Writing all the 
different terms explicitly gives 

a 



B{ki,k2, ks) = n 



M 

— I (7(fcl)[/(fc2)(/(fc3) 
P 



+ ^yU.(fcl)^^) [/(fc2)(7(fc3)f/c(fci) 



+ X]^t/(fcl) f y ) U,{k2)u,{h,)Ucik-,) 



Ms 



p2 ■ ~ ■ p 



U{ki,)Uc{ki)Uc{k2) 



+EEE 



miirijmk 



Ui{ki)uj{k2)uk{k'i) 



X f/e(fcl)[/c(fe2)C/c(fc3) 



+y y — «,(fci) ( ^ ) u,(k2)uj{hi)u,{k^) 



E 



Ui(k-i_)ui{k2)ui{k-i) 



(27) 



where kz ~ — fei — k2. Further integrations over appropri- 
ate window functions yield the spatially smoothed skewness 
(see, e.g., Scoccimarro et al. 2001; Takada & Jain 2002; 
Cooray & Sheth 2002). Similar relations hold for the higher 
order correlations. 

In the same small mass and size limits which we used 
when studying ^(r) and P{k), the bispectrum is dominated 
by the first and last terms of equation (p7|). 



4 AVERAGING OVER NUMBERS AND 
MASSES 

So far we have assumed that the masses and numbers of sub- 
clumps were fixed. This section shows the result of allowing 
for a distribution of masses and numbers. Expressions for 
the power spectra are much simpler than for (,{r), so we will 



only consider P{k) here. Analogously, the averaged bispec- 
trum B is much simpler to evaluate than ^123. For what is to 
follow, it is useful to write explicitly that the quantity in the 
previous section is computed at fixed values of M, M and 
m = {mi, • • • , rriAf}, where the masses and density profiles 
of the N subclumps may be different. The quantity we are 
after is 



dM 



dN(M) 
dM 



p{N\M)p{m\N, M) P{k\m,N, M),(28) 



where dN/ dM is the number density of parent M halos (usu- 
ally called the universal halo mass function), p[N\M) is the 
probability an M halo has N subclumps, p{m\N , M) is the 
probability that the N subclumps were {mi, • • • , mj^f} given 
that there were N of them in the M halo, P{k\m,M , M) is 
the quantity we computed in the previous section, and the 
sum is over all values of A/" and m. To proceed, we need 
models of these different distributions. 

If all the subclumps are identical, and there are A/" of 
them, then each of the terms in the sums over i and j are the 
same, so the contribution to the power from the smooth- 
subclump cross-terms is A/" times the contribution from 
a single clump, and the contribution from the subclump- 
subclump terms is oc Af{Af - l)Ul. In this case, the total 
power depends on the first two factorial moments of the 
p(N\M) distribution. If the subclumps are much smaller 
than their parents, then equation (M) implies that 



(29) 



P{k) ^ ^(y) |l-A/'('n/M)"| U(kf 



+ AAn ( — ) u{ky 



The factor Nn is simply the number density of the sub- 
clumps, so it may be useful to think of the two terms above 
separately. The first term requires knowledge of the mass 
function of the parent halos, whereas the second requires 
the 'mass function' of the subclumps. As discussed in Sec- 
tion ^, accurate formulae for the parent mass function are 
available; subclump mass functions are only just becoming 
available. 

In practice, halos of the same M have different sub- 
structure distributions. For example, we expect that there 
will be some scatter around the mean number of subclumps 
Nci = {N\M) in an M— halo, as well as in the fraction of 
mass associated with the clumped component: Md/M. How- 
ever, if the scatter around the typical subclump configura- 
tion is small, then we should get a reasonable estimate of the 
effects of substructure if we use a good model of the typi- 
cal configuration, and neglect the fact that there is actually 
some scatter around it. This is our strategy. 

Let dn{m\M)/dm denote the typical subclump mass 
function. Then 



P{k) = j , 



+ 2 



X [ dm^^^^^^^u{k\m)UM 
J dm p 
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/ 



dM 



dM J dm\ p 

dn(m2|M, mi) m2 

XM(fc|m2, M, mi) Uc{k\mi,m2) 
dN(M) 
dM 



- u{k\mi, M) 



, dn(m\M) m ,,. ,^,2 /o^\ 
xdm — ^-i — - u(k\m, M) . (30) 
dm p 

The third term is the tricky one: it requires the joint distri- 
bution of mi- and m2-subclumps in an M— halo. Following 
Sheth & Lemson (1999), we have chosen to write this as 
the typical number of mi-subclumps in an M— halo times 
the typical number of m2-subclumps in M-halos which have 
an mi-subclump. Similarly, the notation f7c(fe|mi, m2) is in- 
tended to allow for the possibility that the distribution of 
m2-subclumps within the parent M-halo may depend on 
whether or not the parent contains an mi-subclump. 

In the limit in which the subclumps are much smaller 
than the parent, and contribute a small fraction of the total 
mass, the expression above becomes 



P{k) 



I 



dN{M) 
dM 

dN{M) 



U{k\Mf 

, dn(m\M) m^ , , ,,2 
dm — V — - u(k\m, MY. 
dm P 



The integral over m. is from zero to M. If u(fc|m, M) is only 
a function of m (i.e., if the density profile of an m subclump 
is independent of the mass of the parent in which it sits), 
then the order of the integrals above can be rearranged. If 
wo define the mass function of the subclumps as 



diVc(m) 
dm 

then 

P(k) ^ 



f 

J m 



dM 



dN(M) dn{m\M) 



dM 



dm 



(31) 



(32) 



In this limit, the total power is well approximated by adding 

to the power associated with the smooth parent halos the 
contribution which comes integrating over the subclump 
mass function. 



5 DETAILS 

The results presented above are general. When used to 
model large-scale structure simulations, one will almost al- 
ways use models of the halo density profiles and halo mass 
functions which have been found to provide good descrip- 
tions of numerical simulations. Specifically, the density pro- 
files F and / are usually described using the functional form 
given by Navarro, Prenk & White (1997), 



p{r 



Avir 



3Q (r/r,)(l -h r/r^)2 



(33) 



where the profile is truncated at the virial radius rvir, 
c = rvir/'f's is called the concentration, and /(c) = l/[ln(l -|- 
c) — c/(l -|- c)]. The virial radius Tvir is defined by requir- 
ing that m = 47rr^ijP Avir. For spatially fiat universes with 
^0 = (1, 0.3) and A = 1 - f^, Avir = (178, 340). The Fourier 
transform of the density run around such a halo of mass m 
is 



u{k\m) = f{c) 



lK(si[K(l-|-c)] -si(K; 



+ COS «:( Ci[K(l -I- c)] - C\{n) 



sin(Kc) 

k(i + c; 



(34) 



(Scoccimarro et al. 2001), where n = krvii/c, Si{x) = 

dt sm{t)/t is the sine integral and Ci(a;) = — dt cos{t)/t 
is the cosine integral function. The concentration parameter 
of the halos depends on halo mass; we use the parametriza- 
tion of this dependence given by Bullock et al. (2001): 



c(m) 



1 + z 



4^) 



(35) 



Tha parent halo mass function is well described by 



dN(M, z) dM ^ , , dv 
p dM M ■'^ ' v' 



(36) 



where v ■ 
and 



(5gc(2)/c^(Af), p is the background mass density. 



ufiu) = A(p) (1 + {quy^) (g) exp (-|^) , (37) 

with p w 0.3, A{p) = [1 + 2-*'r(l/2 -p)/V7r]"^ w 0.3222, 
and q « 0.75 (Sheth & Tormen 1999). Here sj{z) is the crit- 
ical density required for spherical collapse at z, extrapolated 
to the present time using linear theory, and 



\M) 



4-K 
(2^ 



dk 

^k^Punik) W^ikRo), 



(38) 



where W(a;) = (3/x^)[sin(a:)-a:cos(a:)] and i?o = (3M/47rp)^/^. 
That is to say, cr(M) is the rms value of the initial fluctuation 
field when it is smoothed with a tophat filter of comoving 
size -Ro, extrapolated using linear theory to the present time. 
If p = 1/2 and q = 1, then dN/dM is the same as the univer- 
sal mass function first written down by Press & Schechter 
(1974). Note that the mass function is normalized so that 



dM MdN{M)/dM = p. 



(39) 



The subclump distribution is less well constrained. There- 
fore, we will consider two models. The first is motivated by 
the results of numerical simulations (Tormen et al. 1998; 
Ghigna et al. 1999) which suggest that, on average, a power 
law in mass is a reasonable model of the typical subclump 
mass function: 



dn{m\M) 
dm 



dm ■■ 



/■M\'^ dm 

No[ — ) 

\ m / m 



where A* < 1, (40) 



and A^o is a normalization constant which is set by the frac- 
tion of mass in subclumps: 



M 



f 

Jo 



dm 



m dn{m\M) No 



M dm 



1-M 

In this model, although the total number of subclumps may 
diverge, the mass contained in them does not. The mean 
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square subclump mass, which is related to the total number 
of subclump pairs is 



m\2 dn(m\M) ^ ^2 No 
Q X Af / dm 2 — n 2 ~ fi 

Simulations suggest jj, ~ 0.9 and Md/M ~ 0.1. If Puu(k) oc 
fc" and we use the Press-Schechter form for dN/dM, then in 
this power law subclump model, the subclump mass function 
is 

Tn^ dN^{m) _ 2(''-i)/°' /A/.X^'^ 
~p dm ^ r(l/2) \~) 

-(|-^.^).("> 

where a = (n + 3)/3 (the value n ~ —3/2 is a reason- 
able approximation to most CDM models). At small masses, 
mIMt <^ 1, the subclump mass function is a power law: 
dNc/dm cx m~'^~'^, whereas at m/Af* ^ 1 it drops expo- 
nentially. 

In practice, the simulations only have finite mass resolu- 
tion, so a power law model with a minimum mass cutoff may 
be more appropriate. A low-mass cutoff may also be a use- 
ful model if one is only interested in subclumps which could 
host galaxies. In this case, the mass fraction in subclumps 
is 



M 



m dn(m\M) ,^ 1 
dm — — V-i — = iVo — 
M dm 



1-M 



The lower mass cutoff means that the number of subclumps 
no longer diverges: 



dm 



Therefore, when Mmin <^ M, the first term is the dominant 
one, and the number of subclumps increases as a power law 
in M. In practice, p is sufficiently close to unity that the 
mass fraction in subclumps is approximately No/{l — n) only 
for the most massive halos. 

Our second model of dn{m\AI)/dm is to assume that 
it is given by the progenitor distribution at some higher 
redshift zi , and all progenitors lost the same fraction of mass 
as they fell into the parent they now occupy, so they are all 
a fraction / of the mass they were at zi. In this case. 



dN,{m,zo) ^ dN{M,zo) dn{m/ f, zi\M, zo) 



dm 



dM 



dm 



(42) 



Mass conservation means that this integral reduces to the 
mass function at the higher redshift: dNc{m,zo)/dm = 
dN{m/ f , Zi) / dm. The total amount of mass in these sub- 
clumps is 



dmmdNc{m)/dm — f / dm (m/ f) dN{m/ f) /dm — f p. 



This provides a convenient way of thinking about the sub- 
clump mass function. In practice, it is likely that the sub- 
clumps which have survived to the present come from a dis- 
tribution of earlier epochs, rather than a single epoch zi, 
and the distribution of zi might peak at different redshifts 
for different values of m/M. Ignoring this subtlety, which 
is what the model above does, should be a reasonable ap- 
proximation if the distribution around the mean zi is not 



broad. In practice, for small values of m/M, the subclump 
mass function is reasonably well approximated by a power 
law, and so this model also gives a subclump mass function 
which is in reasonable agreement with the simulations. 

The virtue of this model is that it allows a simple esti- 
mate of the shape of the correlation function in the limit in 
which the subclumps are each a small fraction of the mass of 
the parent halo, and the total mass in subclumps is a small 
fraction of the total mass. Recall that, in this limit, the total 
correlation function is well approximated by the sum of two 
terms. The first term is the contribution from pairs which are 
from the smoothly distributed component, and the second is 
from the subclumps. The subclump contribution, then, can 
be determined by rescaling the contribution from the smooth 
component at zi. If the particles which were stripped from 
the halos present at zi were random particles, then the num- 
ber of pairs is lower by a factor of since the halos at zi 
were a factor of (l + zi)^ denser, the contribution to the cor- 
relation function is different by a factor of /'^(l + zi)'^. This 
yields a factor which is probably slightly smaller than unity. 
On the other hand, if the mass was stripped from the earlier 
halos in shells, much like layers off an onion, as simulations 
suggest, then a better model of the subclump contribution is 
to truncate the halos which were present at zi at a fraction 
~ / of their virial radii when estimating how the number 
of pairs changes with scale. Although this changes the ac- 
tual shape of the correlation function (e.g., the subclump 
pairs are shifted to scales which are a factor of / smaller), 
the typical factor by which the correlations are affected is 
/"^(l^-^l)^ which can be considerably larger than unity. 



6 DISCUSSION 

We have shown how to incorporate the effects of substruc- 
ture into the halo model description of the nonlinear density 
field. Accounting for this substructure is important on scales 
smaller than the virial radii of typical halos. The effects are 
more pronounced for statistics which treat the subclumps 
preferentially, such as the power spectrum measured in stud- 
ies of weak galaxy-galaxy gravitational lensing. Substruc- 
ture will also change the dynamics within halos. Although 
we have not done so here, it is straightforward to insert our 
model for substructure into the halo model of the cosmic 
virial theorem, and the mean pairwise velocity and velocity 
dispersion developed in Sheth et al. (2001). 

The stable clustering limit is a physically appealing de- 
scription of clustering on small scales (Peebles 1980). It has 
been argued that a model with smooth halos is inconsistent 
with this limit (Ma & Fry 2000; Scoccimarro et al. 2001). 
Substructure changes the shape of the small scale power 
spectrum (c.f.. Fig. ^); at least in principle, it can bring 
the halo model predictions into agreement with the stable 
clustering solution. However, it is not obvious that stable 
clustering is, indeed, the correct physical limit. Smith et al. 
(2002) argue that the stable clustering assumption is incon- 
sistent with the results of high resolution numerical simu- 
lations. They also find that the simulations do not follow 
the small scale scaling predicted by models in which halos 
are smooth. Once an accurate model of the subclump mass 
function is available, it will be interesting to compare the 
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predictions of our description of substructure with their re- Tormen G., Diaferio A., Syer D., 1998, MNRAS, 299, 728 

suits. 

Although we have focussed primarily on the implica- 
tions of substructure for the halo model of nonlinear clus- 
tering, our results have a wide range of other applications. 
For example, excess power in the Fourier transforms of im- 
ages of galaxies or distant clusters can be used to infer the 
existence of spiral arms or substructure. This is the sub- 
ject of work in progress. Closely related is the question of 
what images of high redshift galaxies may look like. Obser- 
vations through a filter which has a fixed wavelength range 
probe the emission from high redshift galaxies at shorter 
restframe wavelengths than for galaxies at low redshift. If 
obscuration by dust is not a problem, and the UV lumi- 
nosity is dominated by patchy star forming regions, then 
the images of high redshift galaxies should show consider- 
able substructure. Our results suggest that, in this case, the 
power spectrum obtained by Fourier transforming the image 
of a high redshift patch of sky should show an increase in 
small scale power. 

In addition, although we have phrased the entire dis- 
cussion of substrucutre in terms of spatial statistics, this is 
not really necessary. Large databases describing various ob- 
served characteristics of galaxies are now becoming available 
(e.g., the 2dFGRS and SDSS surveys). If some of n observ- 
ables are correlated with others, the data will not fill the 
full n— dimensional space available: the data set itself can 
be thought of as being clumpy, and the various clumps in 
dataspace may themselves have substructure. The formalism 
developed here provides a way of discovering, quantifying 
and modelling such substructure. 

We thank Masahiro Takada for helpful discussions. B.J. ac- 
knowledges financial support from a NASA-LTSA grant and 
a Keck foundation grant. 
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